import numpy as np

# CCD
BAND_NUV= 0
BAND_u  = 1
BAND_g  = 2
BAND_r  = 3
BAND_i  = 4
BAND_z  = 5
BAND_Y  = 6

CCD_bandnames = ['NUV', 'u', 'g', 'r', 'i', 'z', 'y']

def init_ccd_bands():
    """
    初始化CCD在焦面上对应的波段信息
    """
    BANDs = [BAND_NUV, BAND_NUV, BAND_NUV, BAND_NUV, \
             BAND_u, BAND_u, \
             BAND_g, BAND_g, \
             BAND_r, BAND_r, \
             BAND_i, BAND_i, \
             BAND_z, BAND_z, \
             BAND_Y, BAND_Y, BAND_Y, BAND_Y ]
    return BANDs

def init_ccd_pos_in_focus():
    """
    初始化CCD在焦面上的坐标
    返回值：
        ccd_pos_in_focus
    备注：目前仅考虑成像部分，无缝光谱部分暂未考虑进来
    """
    ccd_pos_in_focus = np.zeros((18,2))
    #----------------NUV  4----------------------
    ccd_pos_in_focus[0][0] =   0.0030/28.
    ccd_pos_in_focus[0][1] =  -0.1475/28.
    ccd_pos_in_focus[1][0] =   0.0030/28.
    ccd_pos_in_focus[1][1] =  -0.0462/28.
    ccd_pos_in_focus[2][0] =  -0.0952/28.
    ccd_pos_in_focus[2][1] =  -0.0462/28.
    ccd_pos_in_focus[3][0] =  -0.0952/28.
    ccd_pos_in_focus[3][1] =   0.0551/28.

    #--------------------u  2---------------------
    ccd_pos_in_focus[4][0] = -0.0952/28.
    ccd_pos_in_focus[4][1] = -0.1475/28.

    ccd_pos_in_focus[5][0] =  0.0030/28.
    ccd_pos_in_focus[5][1] =  0.0551/28.

    #--------------------g  2---------------------
    ccd_pos_in_focus[6][0] =  0.1012/28.
    ccd_pos_in_focus[6][1] = -0.0462/28.

    ccd_pos_in_focus[7][0] = -0.1934/28.
    ccd_pos_in_focus[7][1] = -0.0462/28.

    #--------------------r  2---------------------
    ccd_pos_in_focus[8][0] = -0.1934/28.
    ccd_pos_in_focus[8][1] = -0.1475/28.

    ccd_pos_in_focus[9][0] =  0.1012/28.
    ccd_pos_in_focus[9][1] =  0.0551/28.

    #--------------------i  2---------------------
    ccd_pos_in_focus[10][0] =  0.1012/28.
    ccd_pos_in_focus[10][1] = -0.1475/28.

    ccd_pos_in_focus[11][0] = -0.1934/28.
    ccd_pos_in_focus[11][1] =  0.0551/28.

    #--------------------z  2---------------------
    ccd_pos_in_focus[12][0] = -0.0952/28.
    ccd_pos_in_focus[12][1] = -0.2488/28.

    ccd_pos_in_focus[13][0] =  0.0030/28.
    ccd_pos_in_focus[13][1] =  0.1564/28.

    #--------------------Y  4---------------------
    ccd_pos_in_focus[14][0] =  0.1012/28.
    ccd_pos_in_focus[14][1] = -0.2488/28.

    ccd_pos_in_focus[15][0] =  0.0030/28.
    ccd_pos_in_focus[15][1] = -0.2488/28.

    ccd_pos_in_focus[16][0] = -0.0952/28.
    ccd_pos_in_focus[16][1] =  0.1564/28.

    ccd_pos_in_focus[17][0] = -0.1934/28.
    ccd_pos_in_focus[17][1] =  0.1564/28.

    return ccd_pos_in_focus


def rotateByAxisY(p, theta):
    """
    rotate vector p about Y axis by an angle of theta
    """
    pr = np.zeros(3)
    pr[0] = p[0]*np.cos(theta) - p[2]*np.sin(theta)
    pr[1] = p[1]
    pr[2] = p[0]*np.sin(theta) + p[2]*np.cos(theta)
    return pr

def rotateByAxisZ(p, theta):
    """
    rotate vector p about Z axis by an angle of theta
    """
    pr = np.zeros(3)
    pr[0] = p[0]*np.cos(theta) - p[1]*np.sin(theta)
    pr[1] = p[0]*np.sin(theta) + p[1]*np.cos(theta)
    pr[2] = p[2]
    return pr

def isInRect(pix_vec, p1, p2, p3, p4):
    """
    判断某个指向是否在p1,p2,p3,p4四个点组成的小矩形范围内
    返回值： True or False
    """
    vec = np.zeros((4,3))
    
    # p1
    a1 = p1[0] - pix_vec[0]
    a2 = p1[1] - pix_vec[1]
    a3 = p1[2] - pix_vec[2]

    b1 = p2[0] - pix_vec[0]
    b2 = p2[1] - pix_vec[1]
    b3 = p2[2] - pix_vec[2]

    vec[0][0] = a2*b3 - a3*b2
    vec[0][1] = a3*b1 - a1*b3
    vec[0][2] = a1*b2 - a2*b1

    # p2
    a1 = p2[0] - pix_vec[0]
    a2 = p2[1] - pix_vec[1]
    a3 = p2[2] - pix_vec[2]

    b1 = p3[0] - pix_vec[0]
    b2 = p3[1] - pix_vec[1]
    b3 = p3[2] - pix_vec[2]

    vec[1][0] = a2*b3 - a3*b2
    vec[1][1] = a3*b1 - a1*b3
    vec[1][2] = a1*b2 - a2*b1

    # p3
    a1 = p3[0] - pix_vec[0]
    a2 = p3[1] - pix_vec[1]
    a3 = p3[2] - pix_vec[2]

    b1 = p4[0] - pix_vec[0]
    b2 = p4[1] - pix_vec[1]
    b3 = p4[2] - pix_vec[2]

    vec[2][0] = a2*b3 - a3*b2
    vec[2][1] = a3*b1 - a1*b3
    vec[2][2] = a1*b2 - a2*b1

    # p4
    a1 = p4[0] - pix_vec[0]
    a2 = p4[1] - pix_vec[1]
    a3 = p4[2] - pix_vec[2]

    b1 = p1[0] - pix_vec[0]
    b2 = p1[1] - pix_vec[1]
    b3 = p1[2] - pix_vec[2]

    vec[3][0] = a2*b3 - a3*b2
    vec[3][1] = a3*b1 - a1*b3
    vec[3][2] = a1*b2 - a2*b1

    flag = True
    scale = 1e8

    flag1 = ( vec[0][0]*vec[1][0] 
            + vec[0][1]*vec[1][1] 
            + vec[0][2]*vec[1][2] ) * scale
    for i in range(1,4):
        flagn = ( vec[i][0]*vec[(i+1)%4][0] 
                + vec[i][1]*vec[(i+1)%4][1] 
                + vec[i][2]*vec[(i+1)%4][2] )*scale
        if flag1 * flagn < 0:
            flag = False
            break

    return flag

def computeCCDPos(lon, lat, y, z):
    """
    计算单个CCD四个角在天球上的指向（3D坐标）
    """
    ccd_y = 0.00324295  # Unit: radian
    ccd_z = 0.00314273  # Unit: radian

    py, pz = y, z
    p1 = np.zeros(3)
    p1[0] = np.sqrt( 1 - py*py - pz*pz )
    p1[1] = py
    p1[2] = pz
    p11 = rotateByAxisY( p1, np.radians(lat) )
    P1  = rotateByAxisZ( p11, np.radians(lon) )

    py = y + ccd_y
    pz = z
    p2 = np.zeros(3)
    p2[0] = np.sqrt( 1 - py*py - pz*pz )
    p2[1] = py
    p2[2] = pz
    p22 = rotateByAxisY( p2, np.radians(lat) )
    P2 = rotateByAxisZ( p22, np.radians(lon) )

    py = y + ccd_y
    pz = z + ccd_z
    p3 = np.zeros(3)
    p3[0] = np.sqrt( 1 - py*py - pz*pz )
    p3[1] = py
    p3[2] = pz
    p33 = rotateByAxisY( p3, np.radians(lat) )
    P3 = rotateByAxisZ( p33, np.radians(lon) )

    py = y
    pz = z + ccd_z
    p4 = np.zeros(3)
    p4[0] = np.sqrt( 1 - py*py - pz*pz )
    p4[1] = py
    p4[2] = pz
    p44 = rotateByAxisY( p4, np.radians(lat) )
    P4 = rotateByAxisZ( p44, np.radians(lon) )
    return P1, P2, P3, P4
